#THIS SCRIPT SETS UP THE DATA FOR ANALYSIS, BY CREATING SIMULATIONS OF IDEAL POINTS
rm(list=ls(all=TRUE))
setwd('~/Dropbox/Move_the_median/Replication Data')

library(foreign)
library(arm)
library(plyr)
library(car)
library(reshape2)
options(digits=2)
options(scipen = 100)

n.sims <- 1000
#save number of sims for other script
save.image("n.sims.image")

#read in data
nom.data <- read.dta("nominee_level_data_updated_for_replication.dta", convert.underscore=T)
nom.data <- nom.data[order(nom.data$date.vote),]
head(nom.data)

rc.data <- read.dta("roll_call_data_updated_for_replication.dta", convert.underscore=T)
rc.data <- rc.data[order(rc.data$year.announced),]
pres.codes <- read.csv("pres_party_by_congress.csv") #for filibuster pivot
pres.codes <- subset(pres.codes, cong>=70)

#######################
#now, sequentially set up make matrix of n.sim simulations for each actor's ideal points, for each nominee
###################################################################################################################
#1) Senate median
###################################################################################################################
	#DW 
sen.data <- read.dta("NOMINATE_DW_from_senate_estimates.dta")
sen.data <- subset(sen.data,  statenm!="USA" & cong >=70) #Drop presidents
sen.data <- subset(sen.data, select=c(cong, idno, dwnom1, bootse1))

data.temp <- cbind(sen.data,t(mapply(rnorm,n.sims,sen.data$dwnom1,sen.data$bootse1)))#for each member of congress,
                                                                                     #simulate their ideal points 1,000 times 
                                                                                      #based on their ideal point and standard error
data.temp$dwnom1 <- NULL
data.temp$bootse1 <- NULL
data.reshaped <- melt(data.temp, id.vars=c("cong", "idno"), #reshape such that unit is member/sim
                      variable.name="simulation", value.name="simulated.DW") 

#now take median in each cong, for each simulation
un.cong <- unique(data.reshaped$cong)
sen.median.matrix <- matrix(NA,nrow=length(un.cong), ncol=n.sims)

for (i in 1:length(un.cong)){
 	ok <-   data.reshaped$cong==un.cong[i] 
 	sen.median.matrix[i,] <- tapply(data.reshaped$simulated.DW[ok], data.reshaped$simulation[ok], median )#get median in each sim
	#colnames(sen.median.matrix[i,] ) <- i
}

sen.median.matrix <- data.frame(sen.median.matrix)
rownames(sen.median.matrix) <- un.cong
colnames(sen.median.matrix) <- 1:n.sims
sen.median.matrix$congress <- un.cong
#reshape so that each congress has n.sim observations
cong.reshaped.dw <- melt(sen.median.matrix,  id.vars=c("congress"), variable.name="simulation", value.name="sen.median.sim.dw")
# #now merge onto nominee data 
df <- subset(nom.data, select=c(nominee, congress))
cong.reshaped.dw <- join(df, cong.reshaped.dw)

#FILIBUSTER PIVOT
#NEED TO MERGE PRESIDENT TO FIGURE OUT FILIBUSTER PIVOT
dim(data.reshaped)
data.reshaped <- join(data.reshaped,pres.codes, match="first")
dim(data.reshaped)
head(data.reshaped)

#now take FP in each cong, for each simulation
	#pivot was 2/3 up until and not including the 94th Congress (only nominee during that Congress was Stevens-he was confirmed after rule change)
	 #need to account for pres
un.cong <- unique(data.reshaped$cong)
sen.pivot.matrix <- matrix(NA,nrow=length(un.cong), ncol=n.sims)

for (i in 1:length(un.cong)){
 	ok <-   data.reshaped$cong==un.cong[i] 
	dem.pres <- rep(unique(ifelse(data.reshaped$pres.party[ok]==100,1,0)),n.sims) #make same length of sims for ifelse statements below
	pre.1975 <- rep(unique(ifelse(data.reshaped$cong[ok] < 94,1,0 )), n.sims) #make same length of sims for ifelse statements below
	foo.33 <- tapply(data.reshaped$simulated.DW[ok], data.reshaped$simulation[ok], quantile, probs=c(.3333))
	foo.4 <- tapply(data.reshaped$simulated.DW[ok], data.reshaped$simulation[ok], quantile, probs=c(.4))
	foo.6 <- tapply(data.reshaped$simulated.DW[ok], data.reshaped$simulation[ok], quantile, probs=c(.6))
	foo.67 <- tapply(data.reshaped$simulated.DW[ok], data.reshaped$simulation[ok], quantile, probs=c(.67777))
 	sen.pivot.matrix[i,] <-	ifelse(dem.pres == 1 & pre.1975 ==1, foo.67,
													ifelse(dem.pres == 1 & pre.1975 ==0, foo.6,
													ifelse(dem.pres == 0 & pre.1975 ==1, foo.33,
													ifelse(dem.pres == 0 & pre.1975 ==0, foo.4))))
}

sen.pivot.matrix <- data.frame(sen.pivot.matrix)
rownames(sen.pivot.matrix) <- un.cong
colnames(sen.pivot.matrix) <- 1:n.sims
sen.pivot.matrix$congress <- un.cong
#reshape so that each congress has n.sim observations
temp.pivot.reshaped.dw <- melt(sen.pivot.matrix,  id.vars=c("congress"), variable.name="simulation", value.name="sen.pivot.sim.dw")
# #now merge onto cong.reshaped.dw [created above]
dim(cong.reshaped.dw )
cong.reshaped.dw <- join(cong.reshaped.dw, temp.pivot.reshaped.dw )
dim(cong.reshaped.dw )

##################################################################################
#BAILEY
##################################################################################
sen.data <- read.dta("bailey_senators_updated.dta")
sen.data <- subset(sen.data, select=c(congress, senator_id, senator_bly, senator_bly_se))
data.temp <- cbind(sen.data,t(mapply(rnorm,n.sims,sen.data$senator_bly, sen.data$senator_bly_se)))

data.temp$senator_bly <- NULL
data.temp$senator_bly_se <- NULL
data.reshaped <- melt(data.temp, id.vars=c("congress", "senator_id"), variable.name="simulation", value.name="simulated.bly")

#now take median in each cong, for each simulation
un.cong <- unique(data.reshaped$congress)
sen.median.matrix <- matrix(NA,nrow=length(un.cong), ncol=n.sims)

for (i in 1:length(un.cong)){
 	ok <-   data.reshaped$congress==un.cong[i] 
 	sen.median.matrix[i,] <- tapply(data.reshaped$simulated.bly[ok], data.reshaped$simulation[ok], median )#get median in each sim
}

sen.median.matrix <- data.frame(sen.median.matrix)
rownames(sen.median.matrix) <- un.cong
colnames(sen.median.matrix) <- 1:n.sims
sen.median.matrix$congress <- un.cong

#reshape so that each congress has n.sim observations
cong.reshaped.bly <- melt(sen.median.matrix,  id.vars=c("congress"), variable.name="simulation", value.name="sen.median.sim.bly")
# #now merge onto nominee data 
df <- subset(nom.data, select=c(nominee, congress))
cong.reshaped.bly <- join(df, cong.reshaped.bly)

#FILIBUSTER PIVOT
dim(data.reshaped)
head(data.reshaped)
names(pres.codes)[names(pres.codes) == 'cong'] <- 'congress'#bailey data uses "congress", not "cong", as var name

data.reshaped <- join(data.reshaped,pres.codes, match="first")#merge pres party
dim(data.reshaped)
head(data.reshaped)

#now take FP in each cong, for each simulation
	#pivot was 2/3 up until and not including the 94th Congress (only nominee during that Congress was Stevens-he was confirmed after rule change)
	 #need to account for pres
un.cong <- unique(data.reshaped$congress)
sen.pivot.matrix <- matrix(NA,nrow=length(un.cong), ncol=n.sims)

for (i in 1:length(un.cong)){
 	ok <-   data.reshaped$congress==un.cong[i] 
	dem.pres <- rep(unique(ifelse(data.reshaped$pres.party[ok]==100,1,0)),n.sims) #make same length of sims for ifelse statements below
	pre.1975 <- rep(unique(ifelse(data.reshaped$congress[ok] < 94,1,0 )), n.sims) #make same length of sims for ifelse statements below
	foo.33 <- tapply(data.reshaped$simulated.bly[ok], data.reshaped$simulation[ok], quantile, probs=c(.3333))
	foo.4 <- tapply(data.reshaped$simulated.bly[ok], data.reshaped$simulation[ok], quantile, probs=c(.4))
	foo.6 <- tapply(data.reshaped$simulated.bly[ok], data.reshaped$simulation[ok], quantile, probs=c(.6))
	foo.67 <- tapply(data.reshaped$simulated.bly[ok], data.reshaped$simulation[ok], quantile, probs=c(.67777))
 	sen.pivot.matrix[i,] <-	ifelse(dem.pres == 1 & pre.1975 ==1, foo.67,
													ifelse(dem.pres == 1 & pre.1975 ==0, foo.6,
													ifelse(dem.pres == 0 & pre.1975 ==1, foo.33,
													ifelse(dem.pres == 0 & pre.1975 ==0, foo.4))))
}

sen.pivot.matrix <- data.frame(sen.pivot.matrix)
rownames(sen.pivot.matrix) <- un.cong
colnames(sen.pivot.matrix) <- 1:n.sims
sen.pivot.matrix$congress <- un.cong

#reshape so that each congress has n.sim observations
temp.pivot.reshaped.bly <- melt(sen.pivot.matrix,  id.vars=c("congress"), variable.name="simulation", value.name="sen.pivot.sim.bly")
# #now merge onto cong.reshaped.bly [created above]
dim(cong.reshaped.bly )
cong.reshaped.bly <- join(cong.reshaped.bly, temp.pivot.reshaped.bly )
dim(cong.reshaped.bly )

#join both DW and BLY together
cong.reshaped <-  join(cong.reshaped.dw, cong.reshaped.bly)

########################################################################################################################################################################################
############################################################################################
#now, using original nom data, generate sims for pres, justices and nominee
temp.data <- subset(nom.data,T)
un.nom <- unique(temp.data$nominee)
############################################################################################
#PRESIDENT
############################################################################################
#DW
pres.matrix <- data.frame(matrix(NA,nrow=length(un.nom), ncol=n.sims))
rownames(pres.matrix) <- un.nom
colnames(pres.matrix) <- 1:n.sims
pres.matrix$nominee <- un.nom

for (i in 1:length(un.nom)){
 	ok <-   pres.matrix$nominee==un.nom[i] 
 	pres.matrix[i,1:n.sims] <- rnorm(n.sims, temp.data$pres.dw[ok], temp.data$pres.dw.se[ok])
}
#reshape so that each pres has n.sim observations
pres.reshaped.dw <- melt(	pres.matrix,  id.vars=c("nominee"), variable.name="simulation", value.name="pres.sim.dw")

#REPEAT FOR BAILEY
pres.matrix <- data.frame(matrix(NA,nrow=length(un.nom), ncol=n.sims))
rownames(pres.matrix) <- un.nom
colnames(pres.matrix) <- 1:n.sims
pres.matrix$nominee <- un.nom

for (i in 1:length(un.nom)){
 	ok <-   pres.matrix$nominee==un.nom[i] 
 	pres.matrix[i,1:n.sims] <- rnorm(n.sims, temp.data$pres.bly[ok], temp.data$pres.bly.se[ok])
}
#reshape so that each pres has n.sim observations
pres.reshaped.bly <- melt(	pres.matrix,  id.vars=c("nominee"), variable.name="simulation", value.name="pres.sim.bly")
#join both DW and BLY together
pres.reshaped <-  join(pres.reshaped.dw, pres.reshaped.bly)

############################################################################################
#NOMINEES
#############################################################################################
#PLAN
#FOR BAILEY, SIMPLY REGRESS 1st-year BLY score on Segal-Cover, then project to all nominees
#For MQ, two steps: 
	#1) TRANSFORM MQ VOTING SCORES INTO DW, USING PRES DW (ALL OR Moraski/Shipan Constrained)
	#2) REGRESS 1st year DW score on Segal Cover, then project to all nominees
	
#1st, create n.sims draws of ideololgy
#MQ
data.temp <- subset(nom.data, select=c(nominee,  nom.segal.cover.11,   nom.MQ.1st, nom.MQ.1st.sd))
data.temp <- cbind(data.temp,t(mapply(rnorm,n.sims,data.temp$nom.MQ.1st,data.temp$nom.MQ.1st.sd)))
data.temp$nom.MQ.1st <- NULL
data.temp$nom.MQ.1st.sd <- NULL
nominee.ideal.points.reshaped.mq <- melt(data.temp, id.vars=c("nominee", "nom.segal.cover.11"), 
	variable.name="simulation", value.name="simulated.1st.MQ")

#Bailey
data.temp <- subset(nom.data, select=c(nominee, nom.bly.1st, nom.bly.1st.sd))
data.temp <- cbind(data.temp,t(mapply(rnorm,n.sims,data.temp$nom.bly.1st,data.temp$nom.bly.1st.sd)))
data.temp$nom.bly.1st <- NULL
data.temp$nom.bly.1st.sd <- NULL
nominee.ideal.points.reshaped.bly <- melt(data.temp, id.vars=c("nominee"), variable.name="simulation", value.name="simulated.1st.bly")
#join together
nominee.ideal.points.reshaped <- join(nominee.ideal.points.reshaped.mq, nominee.ideal.points.reshaped.bly)

#now, run n.sims regression of segal cover on simulated scores
nominee.ideal.points.reshaped$nom.proj.sim.MQ <- NA
nominee.ideal.points.reshaped$nom.proj.sim.bly <- NA
#also:
#Now convert MQ TO Dw
#To avoid endogeneity with theory, use ALL presidents. 
nominee.ideal.points.reshaped <- join(	nominee.ideal.points.reshaped, pres.reshaped)
nominee.ideal.points.reshaped$arctan.pres.sim.dw  <- tan(  (nominee.ideal.points.reshaped$pres.sim.dw*pi)/2)
nominee.ideal.points.reshaped$nom.sim.dw <- NA#distinguish between 1st mapping from MQ to DQ and ...
nominee.ideal.points.reshaped$nom.proj.sim.dw <- NA# projection using Segal cover scores
#to select "constrained" nominees from Moraski and Shipan
# keep.pres <- with(	nominee.ideal.points.reshaped, {ifelse(nominee=="warren"| nominee=="brennan" | nominee=="whittaker" | nominee=="stewart" | nominee=="white" | nominee=="goldberg" 
# 	| nominee=="marshall" | nominee=="burger" | nominee=="blackmum" | nominee=="powell" | nominee=="rehnquist1" | nominee=="rehnquist2" | nominee=="scalia" | nominee=="ginsburg_r" | nominee=="breyer" | nominee=="oconnor",1,0)	})

#save arctan coefs for later (tranforming justices)
arctan.constant.vec <- rep(NA, n.sims)
arctan.beta.vec <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nominee.ideal.points.reshaped$simulation==i
	#BAILEY-- IMPLY REGRESS 1st-year BLY score on Segal-Cover, then project to all nominees
	foo <- lm(simulated.1st.bly~nom.segal.cover.11, data=nominee.ideal.points.reshaped, subset=ok)
	nominee.ideal.points.reshaped$nom.proj.sim.bly[ok] <- coef(foo)[1] + coef(foo)[2]*nominee.ideal.points.reshaped$nom.segal.cover.11[ok]
	#MQ: 1st: TRANSFORM MQ VOTING SCORES INTO DW, USING PRES DW 
	foo <- lm(arctan.pres.sim.dw~simulated.1st.MQ, data=nominee.ideal.points.reshaped, subset=ok)
	nominee.ideal.points.reshaped$nom.sim.dw[ok] <-  2/(pi) * atan(coef(foo)[1] + coef(foo)[2]*nominee.ideal.points.reshaped$simulated.1st.MQ[ok])
 	arctan.constant.vec[i] <- coef(foo)[1] 
	arctan.beta.vec[i] <-  coef(foo)[2]
	#now project using segal cover
	foo <- lm(nom.sim.dw~nom.segal.cover.11,, data=nominee.ideal.points.reshaped, subset=ok)
	nominee.ideal.points.reshaped$nom.proj.sim.dw[ok] <- coef(foo)[1] + coef(foo)[2]*nominee.ideal.points.reshaped$nom.segal.cover.11[ok]
}

#to check how they look-nom versus president
cbind ( tapply(nominee.ideal.points.reshaped$nom.proj.sim.dw, nominee.ideal.points.reshaped$nominee,mean) ,
			tapply(nominee.ideal.points.reshaped$pres.sim.dw, nominee.ideal.points.reshaped$nominee,mean))
cbind(tapply(nominee.ideal.points.reshaped$nom.proj.sim.bly, nominee.ideal.points.reshaped$nominee,mean),
			tapply(nominee.ideal.points.reshaped$pres.sim.bly, nominee.ideal.points.reshaped$nominee,mean))

#############################################################################################
#JUSTICES
#for each sim, need to sim MQ score then convert using arctan
#############################################################################################
#For Bailey, just need to calculate median
#For DW, need to do arctransformation again

#DW
MQ.list.mean <-  paste0("seat",1:9,".MQ")
MQ.list.sd <-  paste0("seat",1:9,".MQ.sd")

#loop over 9 justices on existing natural court
for (j in 1:9){
		#pull out data for relevant seat
		temp.data <- data.frame(nominee=nom.data$nominee, 
				# seat1=nom.data$seat1,	seat2=nom.data$seat2,	seat3=nom.data$seat3,	seat4=nom.data$seat4,	seat5=nom.data$seat5,
				# 	seat6=nom.data$seat6, 	seat7=nom.data$seat7,	seat8=nom.data$seat8,	seat9=nom.data$seat9,
				mean=unlist(list( nom.data[ MQ.list.mean[j]], nom.data[ MQ.list.sd[j]] )[[1]]),
				sd=unlist(list( nom.data[ MQ.list.mean[j]], nom.data[ MQ.list.sd[j]] )[[2]]))
		temp.matrix <- matrix(NA,nrow=length(un.nom), ncol=n.sims)
		temp.matrix <- data.frame(temp.matrix)
		rownames(temp.matrix) <- un.nom
		colnames(temp.matrix) <- 1:n.sims
		temp.matrix$nominee <- un.nom
				for (i in 1:length(un.nom)){
						ok <-   temp.matrix$nominee==un.nom[i]
					MQ.sims <- rnorm(n.sims, 	temp.data[ok,"mean"], temp.data[ok,"sd"])
	 				temp.matrix[i,1:n.sims] <- 2/pi * atan(arctan.constant.vec + arctan.beta.vec*MQ.sims)
		}
		#reshape
		temp.reshaped <- melt(	temp.matrix,  id.vars=c("nominee"), variable.name="simulation", value.name="seat.sim")
		names(temp.reshaped)[names(temp.reshaped)=="seat.sim"] <- paste0("seat",j,".sim.dw")
		#now assign seat variable for merging on names below to get new median
 		temp <- nom.data[names(nom.data)==paste0("seat",j) | names(nom.data)=="nominee"]
		temp.reshaped <- join(temp.reshaped, temp)
		assign(paste0("seat",j,".reshaped"), temp.reshaped) #assign temp matrix to "j"th 
	}

dfs <- list(seat1.reshaped, seat2.reshaped, seat3.reshaped, seat4.reshaped, seat5.reshaped,
	 	seat6.reshaped, seat7.reshaped, seat8.reshaped ,  seat9.reshaped )

all.seats.reshaped.dw <- join_all(dfs)
all.seats.reshaped.dw <- all.seats.reshaped.dw[order(all.seats.reshaped.dw$nominee, all.seats.reshaped.dw$simulation),]
head(all.seats.reshaped.dw)

#NOW, GET J4, J5, and J6-easiest to reshape, then can pick out 4th, 5th, and 6th justices
all.seats.reshaped.dw.long <- melt(all.seats.reshaped.dw, measure.vars=c("seat1.sim.dw", "seat2.sim.dw",
	"seat3.sim.dw", "seat4.sim.dw", "seat5.sim.dw", "seat6.sim.dw", "seat7.sim.dw", "seat8.sim.dw", "seat9.sim.dw")
	 , id.vars=c("nominee","simulation"), value.name="seat.sim.dw")
all.seats.reshaped.dw.long2 <- melt(all.seats.reshaped.dw, measure.vars=c("seat1", "seat2",
	"seat3", "seat4", "seat5", "seat6", "seat7", "seat8", "seat9")
	 , id.vars=c("nominee","simulation"), value.name="seat")
all.seats.reshaped.dw.long2$variable <- NULL
all.seats.reshaped.dw.long$seat <- all.seats.reshaped.dw.long2$seat
all.seats.reshaped.dw.long$simulation <- as.numeric(all.seats.reshaped.dw.long$simulation)
all.seats.reshaped.dw.long2$simulation <- as.numeric(all.seats.reshaped.dw.long$simulation)

#confirm order is the same 
	#add 0 to sims less than 10
temp <- ifelse(all.seats.reshaped.dw.long$simulation<10, paste0(0, all.seats.reshaped.dw.long$simulation),all.seats.reshaped.dw.long$simulation)
all.seats.reshaped.dw.long$nom.sim.dummy <- paste(all.seats.reshaped.dw.long$nominee, temp)
temp <- ifelse(all.seats.reshaped.dw.long2$simulation<10, paste0(0, all.seats.reshaped.dw.long2$simulation),all.seats.reshaped.dw.long2$simulation)
all.seats.reshaped.dw.long2$nom.sim.dummy <- paste(all.seats.reshaped.dw.long2$nominee, temp)
table(all.seats.reshaped.dw.long$nom.sim.dummy==all.seats.reshaped.dw.long2$nom.sim.dummy)#CONFIRMED

#reorder
all.seats.reshaped.dw.long <- all.seats.reshaped.dw.long[order(all.seats.reshaped.dw.long$nom.sim.dummy), ]
all.seats.reshaped.dw.long$ranks <- ave(all.seats.reshaped.dw.long$seat.sim.dw,all.seats.reshaped.dw.long$nom.sim.dummy, FUN=function(x) rank(x))

j4.temp <- all.seats.reshaped.dw.long[all.seats.reshaped.dw.long$ranks==4,]
names(j4.temp)[names(j4.temp)=="seat.sim.dw"] <- "j4old.sim.dw"
j4.temp <- subset(j4.temp, select=c(nominee, simulation, j4old.sim.dw, nom.sim.dummy))

j5.temp <- all.seats.reshaped.dw.long[all.seats.reshaped.dw.long$ranks==5,]
names(j5.temp)[names(j5.temp)=="seat.sim.dw"] <- "j5old.sim.dw"
j5.temp <- subset(j5.temp, select=c(nominee, simulation, j5old.sim.dw, nom.sim.dummy))

j6.temp <- all.seats.reshaped.dw.long[all.seats.reshaped.dw.long$ranks==6,]
names(j6.temp)[names(j6.temp)=="seat.sim.dw"] <- "j6old.sim.dw"
j6.temp <- subset(j6.temp, select=c(nominee, simulation, j6old.sim.dw, nom.sim.dummy))

dfs <- list( j4.temp, j5.temp, j6.temp)
j4.j5.j6.temp.dw <- join_all( dfs)
head(j4.j5.j6.temp.dw )

#now get new median (assuming nominee is confirmed)
#merge all.seats.reshaped.dw.long with nominee ideal points
temp.df <- join(all.seats.reshaped.dw.long, nominee.ideal.points.reshaped)
#merge replacing2 variable onto data
temp2 <- subset(nom.data, select=c(nominee, replacing2))
temp.df <- join(temp.df, temp2)
temp.df <- temp.df[order(temp.df$nom.sim.dummy),]

#now get ideal point of vacating justice--need to account for internal appointments
temp.df <- within(temp.df, {
	temp.seat <- ifelse(nominee=="fortas2" & seat=="fortas", "fortas2",
			 ifelse(nominee=="rehnquist2" & seat=="rehnquist", "rehnquist2",
			 ifelse(nominee=="roberts_j" & seat=="rehnquist", "rehnquist2",	seat)))
	j.exit.sim.dw <- ifelse(temp.seat==replacing2, seat.sim.dw, NA)
})

tapply(temp.df$j.exit.sim.dw, temp.df$nominee, mean, na.rm=T)

head(temp.df[order(temp.df$nom.sim.dummy),],50)
temp.exit.dw <- subset(temp.df, !is.na(j.exit.sim.dw))#save and merge below
temp.exit.dw <- subset(temp.exit.dw , select=c(nom.sim.dummy,j.exit.sim.dw))

#now replace ideal of nominee with ideal of vacating justice
temp.df <- within(temp.df, {
	seat.sim.adjusted.j5 <- ifelse(seat==replacing2, nom.proj.sim.dw, seat.sim.dw)
	ranks <- ave(seat.sim.adjusted.j5 ,nom.sim.dummy, FUN=function(x) rank(x))
})

j5new.temp.dw <- temp.df[temp.df$ranks==5,]
names(j5new.temp.dw)[names(j5new.temp.dw)=="seat.sim.adjusted.j5"] <- "j5new.sim.dw"
j5new.temp.dw <- subset(j5new.temp.dw, select=c(nominee,simulation, nom.sim.dummy, j5new.sim.dw))

########################################################################################################################################################################################
#REPEAT FOR BAILEY JUSTICES
########################################################################################################################################################################################
bly.list.mean <-  paste0("seat",1:9,".bly")
bly.list.sd <-  paste0("seat",1:9,".bly.sd")

#loop over 9 justices on existing natural court
for (j in 1:9){
		#pull out data for relevant seat
		#keep only bailed nominees
		temp.data <- subset(nom.data,bailey.observation==1)
		un.nom <- unique(temp.data$nominee)
		temp.data <- data.frame(nominee=temp.data$nominee, 
				# seat1=nom.data$seat1,	seat2=nom.data$seat2,	seat3=nom.data$seat3,	seat4=nom.data$seat4,	seat5=nom.data$seat5,
				# 	seat6=nom.data$seat6, 	seat7=nom.data$seat7,	seat8=nom.data$seat8,	seat9=nom.data$seat9,
				mean=unlist(list( temp.data[ bly.list.mean[j]], temp.data[ bly.list.sd[j]] )[[1]]),
				sd=unlist(list( temp.data[ bly.list.mean[j]], temp.data[ bly.list.sd[j]] )[[2]]))
		temp.matrix <- matrix(NA,nrow=length(un.nom), ncol=n.sims)
		temp.matrix <- data.frame(temp.matrix)
		rownames(temp.matrix) <- un.nom
		colnames(temp.matrix) <- 1:n.sims
		temp.matrix$nominee <- un.nom
				for (i in 1:length(un.nom)){
					ok <-   temp.matrix$nominee==un.nom[i]
					bly.sims <- rnorm(n.sims, 	temp.data[ok,"mean"], temp.data[ok,"sd"])
	 				temp.matrix[i,1:n.sims] <-	bly.sims 
		}
		#reshape
		temp.reshaped <- melt(	temp.matrix,  id.vars=c("nominee"), variable.name="simulation", value.name="seat.sim")
		names(temp.reshaped)[names(temp.reshaped)=="seat.sim"] <- paste0("seat",j,".sim.bly")
		#now assign seat variable for merging on names below to get new median
 		temp <- nom.data[names(nom.data)==paste0("seat",j) | names(nom.data)=="nominee"]
		temp.reshaped <- join(temp.reshaped, temp)
		assign(paste0("seat",j,".reshaped"), temp.reshaped) #assign temp matrix to "j"th 
	}

dfs <- list(seat1.reshaped, seat2.reshaped, seat3.reshaped, seat4.reshaped, seat5.reshaped,
	 	seat6.reshaped, seat7.reshaped, seat8.reshaped ,  seat9.reshaped )

all.seats.reshaped.bly <- join_all(dfs)
all.seats.reshaped.bly <- all.seats.reshaped.bly[order(all.seats.reshaped.bly$nominee, all.seats.reshaped.bly$simulation),]
head(all.seats.reshaped.bly)

#NOW, GET J4, J5, and J6-easiest to reshape, then can pick out 4th, 5th, and 6th justices
all.seats.reshaped.bly.long <- melt(all.seats.reshaped.bly, measure.vars=c("seat1.sim.bly", "seat2.sim.bly",
	"seat3.sim.bly", "seat4.sim.bly", "seat5.sim.bly", "seat6.sim.bly", "seat7.sim.bly", "seat8.sim.bly", "seat9.sim.bly")
	 , id.vars=c("nominee","simulation"), value.name="seat.sim.bly")

all.seats.reshaped.bly.long2 <- melt(all.seats.reshaped.bly, measure.vars=c("seat1", "seat2",
	"seat3", "seat4", "seat5", "seat6", "seat7", "seat8", "seat9")
	 , id.vars=c("nominee","simulation"), value.name="seat")
all.seats.reshaped.bly.long2$variable <- NULL
all.seats.reshaped.bly.long$seat <- all.seats.reshaped.bly.long2$seat
all.seats.reshaped.bly.long$simulation <- as.numeric(all.seats.reshaped.bly.long$simulation)
all.seats.reshaped.bly.long2$simulation <- as.numeric(all.seats.reshaped.bly.long$simulation)

#confirm order is the same 
	#add 0 to sims less than 10
temp <- ifelse(all.seats.reshaped.bly.long$simulation<10, paste0(0, all.seats.reshaped.bly.long$simulation),all.seats.reshaped.bly.long$simulation)
all.seats.reshaped.bly.long$nom.sim.dummy <- paste(all.seats.reshaped.bly.long$nominee, temp)
temp <- ifelse(all.seats.reshaped.bly.long2$simulation<10, paste0(0, all.seats.reshaped.bly.long2$simulation),all.seats.reshaped.bly.long2$simulation)
all.seats.reshaped.bly.long2$nom.sim.dummy <- paste(all.seats.reshaped.bly.long2$nominee, temp)
table(all.seats.reshaped.bly.long$nom.sim.dummy==all.seats.reshaped.bly.long2$nom.sim.dummy)#CONFIRMED

#reorder
all.seats.reshaped.bly.long <- all.seats.reshaped.bly.long[order(all.seats.reshaped.bly.long$nom.sim.dummy), ]
all.seats.reshaped.bly.long$ranks <- ave(all.seats.reshaped.bly.long$seat.sim.bly,all.seats.reshaped.bly.long$nom.sim.dummy, FUN=function(x) rank(x))

j4.temp <- all.seats.reshaped.bly.long[all.seats.reshaped.bly.long$ranks==4,]
names(j4.temp)[names(j4.temp)=="seat.sim.bly"] <- "j4old.sim.bly"
j4.temp <- subset(j4.temp, select=c(nominee, simulation, j4old.sim.bly, nom.sim.dummy))

j5.temp <- all.seats.reshaped.bly.long[all.seats.reshaped.bly.long$ranks==5,]
names(j5.temp)[names(j5.temp)=="seat.sim.bly"] <- "j5old.sim.bly"
j5.temp <- subset(j5.temp, select=c(nominee, simulation, j5old.sim.bly, nom.sim.dummy))

j6.temp <- all.seats.reshaped.bly.long[all.seats.reshaped.bly.long$ranks==6,]
names(j6.temp)[names(j6.temp)=="seat.sim.bly"] <- "j6old.sim.bly"
j6.temp <- subset(j6.temp, select=c(nominee, simulation, j6old.sim.bly, nom.sim.dummy))

dfs <- list( j4.temp, j5.temp, j6.temp)
j4.j5.j6.temp.bly <- join_all( dfs)
head(j4.j5.j6.temp.bly )

#now get new median (assuming nominee is confirmed)
#merge all.seats.reshaped.bly.long with nominee ideal points
temp.df <- join(all.seats.reshaped.bly.long, nominee.ideal.points.reshaped)
#merge replacing2 variable onto data
temp2 <- subset(nom.data, select=c(nominee, replacing2))
temp.df <- join(temp.df, temp2)
temp.df <- temp.df[order(temp.df$nom.sim.dummy),]

#now get ideal point of vacating justice--need to account for internal appointments
temp.df <- within(temp.df, {
	temp.seat <- ifelse(nominee=="fortas2" & seat=="fortas", "fortas2",
			 ifelse(nominee=="rehnquist2" & seat=="rehnquist", "rehnquist2",
			 ifelse(nominee=="roberts_j" & seat=="rehnquist", "rehnquist2",	seat)))
	j.exit.sim.bly <- ifelse(temp.seat==replacing2, seat.sim.bly, NA)
})

tapply(temp.df$j.exit.sim.bly, temp.df$nominee, mean, na.rm=T)

head(temp.df[order(temp.df$nom.sim.dummy),],50)
temp.exit.bly <- subset(temp.df, !is.na(j.exit.sim.bly))#save and merge below
temp.exit.bly <- subset(temp.exit.bly , select=c(nom.sim.dummy,j.exit.sim.bly))

#now replace ideal of nominee with ideal of vacating justice
temp.df <- within(temp.df, {
	seat.sim.adjusted.j5 <- ifelse(seat==replacing2, nom.proj.sim.bly, seat.sim.bly)
	ranks <- ave(seat.sim.adjusted.j5 ,nom.sim.dummy, FUN=function(x) rank(x))
})

j5new.temp.bly <- temp.df[temp.df$ranks==5,]
names(j5new.temp.bly)[names(j5new.temp.bly)=="seat.sim.adjusted.j5"] <- "j5new.sim.bly"
j5new.temp.bly <- subset(j5new.temp.bly, select=c(nominee,simulation, nom.sim.dummy, j5new.sim.bly))

#merge together

temp.df <- subset(nom.data, select=c( nominee,year.announced, date.vote, confirmed, quality, rcvote, bailey.observation))
dfs <- list(temp.df, nominee.ideal.points.reshaped, cong.reshaped, pres.reshaped, all.seats.reshaped.dw, j4.j5.j6.temp.dw, j5new.temp.dw, temp.exit.dw,
		all.seats.reshaped.bly, j4.j5.j6.temp.bly, j5new.temp.bly, temp.exit.bly)

merged.sims <- join_all(dfs)

#indicator for whether the president is above or below the old median (ABOVE==MORE CONSERVATIVE)
merged.sims$above.dw <- ifelse(merged.sims$pres.sim.dw > merged.sims$j5old.sim.dw,1,0)
merged.sims$above.bly <- ifelse(merged.sims$pres.sim.bly > merged.sims$j5old.sim.bly,1,0)

#recode Bailey nominee projected as missing for bailey observations
merged.sims$nom.proj.sim.bly <- ifelse(merged.sims$bailey.observation==0, NA, merged.sims$nom.proj.sim.bly)

#write data to Stata file
write.dta(merged.sims, "simulated_nom_level_data_from_R_new.dta")

#NOW MERGE ONTO ROLL CALL
rc.cut <- subset(rc.data, select=c(nominee, senator.id, senator.dw, senator.dw.se))
#plan: for each senator/nominee pair, gen 1000 draws of senator's ideal point
#DW
rc.cut <- cbind(rc.cut,t(mapply(rnorm,n.sims,rc.cut$senator.dw,rc.cut$senator.dw.se)))
rc.cut$senator.dw<-NULL
rc.cut$senator.dw.se<-NULL
head(rc.cut)
rc.reshaped.dw <- melt(rc.cut,  id.vars=c("nominee", "senator.id"), variable.name="simulation", value.name="senator.sim.dw")
#Baiely
rc.cut <- subset(rc.data, select=c(nominee, senator.id, senator.bly, senator.bly.se))
rc.cut <- cbind(rc.cut,t(mapply(rnorm,n.sims,rc.cut$senator.bly,rc.cut$senator.bly.se)))
rc.cut$senator.bly<-NULL
rc.cut$senator.bly.se<-NULL
head(rc.cut)
rc.reshaped.bly <- melt(rc.cut,  id.vars=c("nominee", "senator.id"), variable.name="simulation", value.name="senator.sim.bly")

#now merge onto to other sen.level data and nom.data
temp <- subset(rc.data, select=c(name,nominee, vote,vote.adjusted, dem.president, sameprty, senator.id))
dfs <- list(temp, merged.sims, rc.reshaped.dw, rc.reshaped.bly)
rc.reshaped <- join_all(dfs)

#write data to Stata file
write.dta(rc.reshaped, "simulated_sen_level_data_from_R_new.dta")



			   

